|
Карань Анна |
|||
| Главная | О себе | Учеба | ФББ МГУ |
Предсказание генов эукариот
Часть 1: Подготовка чтений
Анализ качества чтений
Сначала в этом задании необходимо сделать контроль качества чтении с помощью программы FastQC. Файл с одноконцевыми чтениями я не привожу, он лежит по ссылке /P/y14/term3/block4/SNP/reads/chr9_1.fastq. Как видно из названия файла я исследовала 9 хромосому человека.
| fastqc chr9_1.fastq |
chr9_1_fastqc.html - выходной файл программы в виде html страницы.
Очистка чтений
| java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr9_1.fastq trimm.fastq SLIDINGWINDOW:10:20 MINLEN:50 TRAILING:3 chr9_1.fastq - входной файл trimm.fastq - выходной файл SLINGWINDOW:10:20 - проходит по прочтениям длинной окном длиной 10 и удаляет правый конец прочтения после окна со средним качеством меньше 20 (если такое окно найдется) MINLEN:50 - удаляет прочтения короче 50 TRAILING:3 - удаляют нуклеотиды ниже качества равного 3-м с конца прочтения |
Теперь нужно провести такой же анализ с помощью FastQC после очистки и сравнить с таковым до.
| fastqc trimm.fastq |
trimm_fastqc.html - выходной файл программы в виде html страницы. Далее нужно сравнить чтения до и после очистки. До очистки было 10701 прочтение, после 10272.
На Рис.1. и Рис.2. изображены качества прочтений до и после очистки. Данный график показывает определенные параметры прочтений. Желтые прямоугольники - интетквартильный интервал, т.е. значения от 25% до 75% всех значений, от него сверху и снизу черные линии, это от 10% до 90% значений, остальные на данной графике не отображаются. Красная линия на каждой интерквартиле - медиана, синия линия , проходящая через весь график - среднее значений. Если сравнить эти 2 графика, видно, что после чистки интервал от 10% до 90% значительно сузился, только нижние значения, т.е. с низким качеством, что как раз является целью очистки (на конце более плохие значения, чем в начале на обоих графиках). Также среднее значение после очистки находится внутри интерквартильных вариантов, и медиана возрасла. Если обобщить, то качества значительно улучшилось.
На Рис.3 и Рис.4. изображены распределения числа прочтений в зависимости от качества. И до, и после очистки наибольшее число прочтений приходилось на качество 38-39, однако до очистки встречались прочтения качества, начиная с 2-х, а после очистки лишь с 25, естественно после 25 количество для любого качество совпадает до и после очистки. Это согласуется с заданными параметрами очистки, все низкокачественные чтения были убраны. Per base sequence content и Per sequence GC content почти не изменились, но и до, и после очистки эти графики являются недостоверными (крестик или восклицательный знак на html странице результатов). После очистки Per base sequence content даже ухудшился, возможно из-за уменьшения числа прочтений, хотя оно было незначительным по сравнению с общим числом.
Модуль Kmer исходит из предположения, что любой небольшой фрагмент последовательности не должен иметь позиционное смещение в его библиотеке. Могут быть биологические причины, почему некоторые Kmers обогащены или обедненны в целом, но эти отклонения должны влиять на все позиции в последовательности одинаково. Таким образом, этот модуль измеряет количество каждого 7-меров на каждой позиции в библиотеке, а затем использует биномиальное распределение для поиска существенных отклонений от равномерного покрытия на всех позициях. Любые Kmers с позиционным смещенением обогащения показываются. На Рис.4. показано распределение TACTGAA, значит, остальные распределены равномерно. После очистки перераспределенных kmer не найдено. Любые индивидульно сверхпредставленные последовательности, даже если не присутствуют в достаточно больших количествах, чтобы вызвать модуль сверхпредставленных последовательностей, будут появляться как "острые шипы" обогащения на графике. Возможно эти пики исчезли после очистки, так как был использован параметр SLINGWINDOW:10:20, и часть последовательсностей была удалена, так что их число уже не достигало сверхраспределенности даже для отдельных пиков.
Часть 2: Картирование чтений
Картированиие чтений
Необходимо откартировать чтения с помощью программы BWA. С помощью следующей команды была проиндексирована референсная последовательность - хромосома 9.
| bwa index chr9.fasta |
Далее было построено выравнивание прочтений и референса в формате .sam. align.sam
| bwa mem chr9.fasta chr9_1.fastq |
Анализ выравниваний
Далее выравнивание чтений с референсом переводится в бинарный формат .bam.
| samtools view -b align.sam -o align.bam -b - параметр, показывающий, что формат выходящего файла - bam. -o - выходящий файл |
Далее отсортируем выравнивание чтений с референсом (получившийся после картирования .bam файл) по координате в референсе начала чтения
| samtools sort align.bam align_sort.bam |
Далее проиндексируем отсортированный .bam файл.
| samtools index align_sort.bam |
Далее посмотрим, сколько чтений откартировалось на геном.
| samtools idxstats align_sort.bam > number |
Выходной файл - number. Откартировалось 10692 из 10701 чтений, т.е. не откартировалось 9 чтений.
Часть 3: Анализ SNP
Поиск SNP и инделей
Сначала был создан файл с полиморфизмами в формате .bcf.
| samtools mpileup -uf chr9.fasta align_sort.bam -o pol.bcf |
Далее был создан файл со списком отличий между референсом и чтениями в формате .vcf.
| bcftools call -cv pol.bcf -o pol.vcf |
Было получено 105 SNP и 6 инделей. Среднее значение глубины покрытия среди полиморфизмов - 13,22857, а среднее значение качества чтений - 83,15374. Максимальное значение покрытия - 98, а качества - 225,544. Покрытие - это то, сколько чтений содержит именно это место. Мне пока трудно сказать, какие значения считаются хорошими, а какие нет, из-за отсутствия опыта. Далее нужно описать три полиморфизма из файла pol.vcf.
| Таблица 1. Информация о полиморфизмах | |||||
| Координата | Тип полиморфизма | В референсе | В чтениях | Глубина покрытия | Качество чтений |
| 4985879 | Замена | A | G | 30 | 225.009 |
| 5122854 | Делеция | agg | ag | 12 | 96.4609 |
| 136131056 | Делеция | CGGG | CGG | 5 | 29.4719 |
Аннотация SNP
С помощью скрипта not_indel.py были убраны строки с инделями и получен файл pol_without_indel.vcf. Далее необходимо получить файл, с которым умеет работать программа annovar, сделать это можно с помощью скрипта convert2annovar.pl следующей командой.
| perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 pol_without_indel.vcf -outfile pol_without_indel.avinput |
pol_without_indel.avinput далее используется для аннотации с помощью баз данных: refgene, dbsnp, 1000 genomes, Gwas, Clinvar. 1) База refgene
| perl /nfs/srv/databases/annovar/annotate_variation.pl -out pol.refgene -build hg19 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb |
Были получены файлы pol.refgene.exonic_variant_function (файл с описанием SNP, попавших в экзоны), pol.refgene.log (файл с описанием прохождения работы), pol.refgene.variant_function (файл с описание SNP по группам). В экзоны попали 15 SNP. (Таблица 2)
| Таблица 2. Информация о SNP в экзонах | |||||||
| Координата | Тип замены | Ген | В референсе | В чтениях | Гетеро или гомозиготная замена | Качество чтений | Глубина покрытия |
| 5081780 | синонимичная | JAK2 | G | A | гомо | 221.999 | 50 |
| 5126443 | несинонимичная | JAK2 | T | A | гетеро | 207.009 | 38 |
| 6253571 | синонимичная | IL33 | С | T | гетеро | 225.009 | 56 |
| 136131188 | неизвестно | неизвестно | С | T | гетеро | 13.2188 | 10 |
| 136131315 | неизвестно | неизвестно | С | G | гетеро | 203.009 | 26 |
| 136131322 | неизвестно | неизвестно | G | T | гетеро | 179.009 | 26 |
| 136131415 | неизвестно | неизвестно | С | T | гетеро | 134.008 | 25 |
| 136131461 | неизвестно | неизвестно | G | A | гетеро | 157.009 | 25 |
| 136131592 | неизвестно | неизвестно | G | C | гетеро | 36.0081 | 7 |
| 136131651 | неизвестно | неизвестно | G | A | гетеро | 4.12848 | 3 |
| 136132873 | неизвестно | неизвестно | T | C | гетеро | 225.009 | 42 |
| 136133506 | неизвестно | неизвестно | A | G | гомо | 221.999 | 29 |
| 136135237 | неизвестно | неизвестно | A | G | гомо | 221.999 | 20 |
| 136135238 | неизвестно | неизвестно | T | C | гомо | 221.999 | 19 |
| 136136770 | неизвестно | неизвестно | A | C | гомо | 222.332 | 17 |
| Таблица 3. Информация о SNP в различных группых | ||||
| В интронах | В экзонах | UTR3 | Гомо | Гетеро |
| 79 | 15 | 9 | 66 | 39 |
Из Таблицы 3 следует, что больше всего SNP в интронах, так как замена в некодирующей последовательности не приведет к замене аминокислоты, миксимум, что может произойти в интронах - это нарушатся консенсусные последовательности для их вырезания во время сплайсинга, но они малы, так что это очень маловероятно. SNP ы экзонах были лишь в генах JAK2 и IL33, куда попали SNP. JAK2 - ген янус киназы 2 (так названа из-за двух киназных доменов в одной молекуле), это тирозин-киназа, фосфорилирующие STAT-факторы. IL33 - ген интерлейкина 33, он экспрессируется многими клетками организма, его уровень строго коррелирует с уровнем воспаления в ткани. В отличие от провоспалительного интерлейкина 1 интерлейкин 33 обладает иммунорегуляторными свойствами. Большое число SNP в интронах было в генах белков, отвечающих за группу крови системы ABO. 2) База dbsnp
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out pol.dbsnp -build hg19 -dbtype snp138 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb |
Были получены файлы pol.dbsnp.hg19_snp138_dropped (SNP c rs), pol.dbsnp.hg19_snp138_filtered (без rs), pol.dbsnp.log. 99 SNP имеют rs, 7 нет. 3) База 1000 genomes
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out pol.1000 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb |
Были получены файлы pol.1000.hg19_ALL.sites.2014_10_dropped, pol.1000.hg19_ALL.sites.2014_10_filtered, pol.1000.log. Состав файлов здесь такой же, как в прошлой базе данных. Однако здесь аннотировано лишь 95, а 11 нет. Среднее значение частоты 0,404849133, минимальное - 0,00319489, максимальное - 1. 4) База Gwas
| perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -build hg19 -out pol.gwas -dbtype gwasCatalog pol_without_indel.avinput /nfs/srv/databases/annovar/humandb |
Были получены файлы pol.gwas.hg19_gwasCatalog, pol.gwas.log. В основном файле содежаться SNP, имеющие какое-то значение с медицине. В моем случае 8 SNP аннотированы в данной базе данных. Один отвечает за предрасположенность к болезни Крона, другой эндометриозу, еще один малярии. Есть 3, отвечающих за уровень коагуляции, также за концентрацию гемоглобина, уровень d-димеров, соответственно венозный тромбоэмболизм. 5) База Clinvar
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out pol.clinvar -dbtype clinvar_20150629 -buildver hg19 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb |
Были получены файлы clinvar.hg19_clinvar_20150629_dropped, clinvar.hg19_clinvar_20150629_filtered, pol.clinvar.log. Аннотированы 3 SNP в генах, связанных с группами крови ABO. Итоговый файл по всем SNP - SNP.xlsx.
| Таблица 4. Итоговая таблица по всем командам | |
| Команда | Функция |
| fastqc chr9_1.fastq | Контроль качества чтений |
| java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr9_1.fastq trimm.fastq SLIDINGWINDOW:10:20 MINLEN:50 TRAILING:3 | Очистка чтений |
| bwa index chr9.fasta | Индексация референсной последовательности |
| bwa mem chr9.fasta chr9_1.fastq | Выравнивание прочтений и референса в формате .sam |
| samtools view -b align.sam -o align.bam | Перевод выравнивания в бинарный формат .bam |
| samtools sort align.bam align_sort.bam | Сортировка выравнивания по координате в референсе начала чтения |
| samtools index align_sort.bam | Индексация отсортированного .bam файла |
| samtools idxstats align_sort.bam > number | Подсчет числа чтений откартированных на геном |
| samtools mpileup -uf chr9.fasta align_sort.bam -o pol.bcf | Создание файла с полиморфизмами в формате .bcf |
| bcftools call -cv pol.bcf -o pol.vcf | Создание файла со списком отличий между референсом и чтениями в формате .vcf |
| perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 pol_without_indel.vcf -outfile pol_without_indel.avinput | Получение файла в форме, пригодной для работе в annovar |
| perl /nfs/srv/databases/annovar/annotate_variation.pl -out pol.refgene -build hg19 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb | Аннотация SNP по базе refgene |
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out pol.dbsnp -build hg19 -dbtype snp138 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb | Аннотация SNP по базе dbsnp |
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out pol.1000 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb | Аннотация по базе 1000 genomes |
| perl /nfs/srv/databases/annovar/annotate_variation.pl -regionanno -build hg19 -out pol.gwas -dbtype gwasCatalog pol_without_indel.avinput /nfs/srv/databases/annovar/humandb | Аннотация по базе Gwas |
| perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out pol.clinvar -dbtype clinvar_20150629 -buildver hg19 pol_without_indel.avinput /nfs/srv/databases/annovar/humandb | Аннотация по базе Clinvar |
©Карань Анна, 2015